Stabilization of bright solitons and vortex solitons in a trapless three-dimensional 
Bose-Einstein condensate by temporal modulation of the scattering length 
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Using variational and numerical solutions of the mean-field Gross-Pitaevskii equation we show 
that a bright sohton can be stabilized in a trapless three-dimensional attractive Bose-Einstein con- 
densate (BEG) by a rapid periodic temporal modulation of scattering length alone by using a 
Feshbach resonance. This scheme also stabilizes a rotating vortex soliton in two dimensions. Apart 
from possible experimental application in BEG, the present study suggests that the spatiotemporal 
solitons of nonlinear optics in three dimensions can also be stabilized in a layered Kerr medium with 
sign-changing nonlinearity along the propagation direction. 

PACS numbers: 03.75.Kk, 03.75.Lm 



I. INTRODUCTION 

Solitons are solutions of wave equation where local- 
ization is obtained due to a nonlinear attractive interac- 
tion. Solitons have been noted in optics 1], high-energy 
physics and water waves |3|, and more recently in Bose- 
Einstein condensates (BEC) 3,4]. The bright solitons of 
attractive BEC represent local maxima jj 13, te] , whereas 
dark solitons of repulsive BEC represent local minima [3 . 

A classic text-book example of soliton appears in 
the following one-dimensional nonlinear free Schrodinger 
equation in dimensionless units [3ilil 






*(y,0=0. 



(1) 



The solitons of this equation are localized solution due 
to the attractive nonlinear interaction — |vl'(y,i)p with 
wave function at time t and position y: '^(y,t) — 

\/2\£\ exp{—i£t)sech{y ^/\£\) , with £ the energy 0. The 
Schrodinger equation with a nonlinear interaction — |^P 
does not sustain a localized solitonic solution in two (2D) 
or three dimensions (3D). However, a radially trapped 
and axially free version of this equation in 3D does sus- 
tain such a bright solitonic solution |^ Q which has been 
observed experimentally in BEC 3 . 

To generate a soliton without trap the repulsive kinetic 
pressure has to balance the attractive force. For a con- 
densate of size L, kinetic energy is proportional to L^^ 
whereas attraction is proportional to L^^ in D dimen- 
sions. The effective potential, which is a sum of these two 
terms, has a confining minimum only for D = 1 leading to 
a stable bound state 0, I3 • Thus no stabilization can be 
obtained in 2D or 3D and any attempt to create a soliton 
either leads to collapse or an expansion to infinity. 

A scheme of stabilization of a soliton in two dimensions 
has been suggested 0, [Q recently by a rapid periodic 
temporal modulation of scattering length a of angular 
frequency w via a —^ ao[l — cshi(ujt)] where uq and c 
are constants and t is time. Such a modulation of the 
scattering length is possible by manipulatin g an external 
magnetic field near a Feshbach resonance |11IJ and has 
been employed in different studies of the BEC U3|- -By 
such a modulation the atomic interaction can be easily 
switched between attractive and repulsive thus resulting 
in a rapid contraction and expansion of the condensate. 
If the constants of modulation are appropriately chosen 
this leads to a stabilization of the condensate in 2D with 



breathing oscillation ^]. However, it has been ''proved" 
by analytic and numerical calculations in Ref. ^Sl that 
such a stabilization does not take place in 3D. 

As the problem of stabilization of a soliton in a trapless 
condensate is of utmost interest in several areas, e.g., op- 
tics 2, 13], nonlinear physics 2] and Bose-Einstein con- 
densate, we revisit this problem and find that a temporal 
modification of the scattering length can also lead to a 
stabilization of the trapless soliton in 3D. We use both 
variational as well as numerical solutions of the mean- 
field time-dependent Gross-Pitaevskii (GP) equation to 
establish our claim. To the best of our knowledge this is 
the first suggestion of stabilization of a trapless soliton 
in 3D. We find numerically that an untrapped attractive 
condensate can maintain a reasonably constant spatial 
profile over a large interval of time through the tempo- 
ral modulation of the s-wave scattering length. We also 
point out a possible reason for the failure to stabilize a 
bright soliton in 3D in Ref. [1(1 ]. 

The present approach is also extended to stabilize vor- 
tex solitons in 2D 14] with angular momentum h per 
atom in the axial (azimuthal) direction. The two dimen- 
sional geometry can be achieved in an axially symmetric 
configuration by applying a strong trap in the axial direc- 
tion Q . This is also equivalent to applying a very weak 
trap in the radial direction. In both cases the radial di- 
mension of the condensate is much larger than the axial 
dimension and a 2D treatment can be justified. Vortex 
solitons are rotating solitons of an attractive condensate 
and it is suggestive that a similar scheme can also be 
used to stabilize a vortex soliton of a trapless conden- 
sate in 3D. After the experimental observation jj^ of 
a vortex in a rotating BEC and theoretical prediction 
J3 of a radially trapped bright vortex soliton, the ex- 
perimental stabilization of trapless vortex solitons seems 
viable. However, the stabilization of a vortex soliton in 
3D calls for a full three-dimensional calculation and is 
beyond the scope of the present study. In the present 
study we stick to a two-dimensional circularly symmetric 
configuration. It is well-known that in a real 3D problem, 
the vortices are unstable against azimuthal perturbation 
which breaks the azimuthal symmetry and calls for a full 
3D treatment. However, these degrees of freedom are ex- 
pected to be partially suppressed in the limit of a very 
strong azimuthal trap or a very weak radial trap when the 
vortex dynamics becomes essentially two dimensional. In 
that limit a circularly symmetric 2D calculation for a 
bright vortex soliton should be sufficient. In this paper 
we find that such a 2D vortex is stable against radial 



perturbation. However, we have not established its sta- 
bility under transverse perturbation. The stability under 
transverse perturbation can be tested by a calculation in 
cartesian coordinates and would be an investigation of 
future interest. 

In Sec. II we present the mean-field model which we 
use in our study. In Sec. Ill and IV, respectively, we 
present the variational and numerical results of our in- 
vestigation and in Sec. V we present our conclusions. 



II. MEAN-FIELD MODEL 

We use the mean-field GP equation for the present in- 
vestigation 13 • In terms of an external reference angular 

frequency Jl and length / = ^Jh|(m^t?) the GP equation 
for the time-dependent Bose-Einstein condensate wave 
function $(r; t) at positio n r and time t can be rewritten 
in dimensionless form as |14lj 









^2 



^1 



8\/27rn I * (f;t) I ^ 



*(f;t) = 0. 



(2) 



where length, time, V^, and wave function are expressed 
in units of //\/2, f^-^, {l|^/2Y'^, and (^/V2)~3/^ respec- 
tively. Here nonlinearity n = Na/l, m is the mass and 
N the number of atoms in the condensate, and a the 
atomic scattering length. The scattering length a and 
nonlinearity n are negative for an attractive condensate 
and positive for a repulsive condensate. In Eq. |(2J) there 
is an axially symmetric harmonic trap with angular fre- 
quency ujp in the radial direction p and lu^ in the axial 
direction z. The normalization condition in Eq. ^ is 
/df|^(f;i)|2 = l. 

The quasi 2D limit of Eq. |(2J is achieved by considering 
Cl ~ ujz » ojp- This condition is satisfied by taking the 
limit Wp — > for a fixed ^ = Uz- This corresponds to a 
pancake-shaped condensate and we look for solution of 
the form ^(r;t) = A{z)il}{p'Tt) with A{z) satisfying the 
one-dimensional ground-state oscillator equation 



1 



(iz2 



A{z) + -rM'A = 7:-4(^)> 



yi7(2^exp[-zV2] 



(3) 



and 



with |^(z)|2 

j^^\A{z)\'^dz = 1. Multiplying Eq. |2Il by A*{z) 
and integrating over z we get the quasi two-dimensional 
GP equation for -tlilp^t): 8] 



d 



-K+ 



4^(2^ 



AnV2TT\^P{p;t)y 



iA(p;t)=0, (4) 



with normalization J dp\ip{p;t)\'^ = 1. In a quantized 
vortex state 2^, with each atom having angular mo- 
mentum Lh along the z axis, ip{p,t) = ip{p,t) exp{iL(J3) 
where ip is the azimuthal angl e. T hen the radial part of 
the GP equation Q becomes p^ 



AnV2T:\ip{p;t)\' 



ip{p;t) = 0, 



(5) 



with normalization 27rJ„ dpp\ip{p]t)\'^ = 1. In Eq. (jsj 
we have introduced a strength parameter d{t) with the 
radial trap. Normally, in the presence of the radial trap 
d{t) = 1. When the radial trap is switched off d{t) will 
be reduced to 0. 

The spherically-symmetric limit of the three- 
dimensional GP equation Q is obtained by taking 
^z = Wp = cUr- In the spherically symmetric configura- 
tion, ^(r,t) — (p{r,t). Then the radial part of the GP 
equation |2j becomes 

_ ■d__^_2d_ u4_ 2 
*a-t a„2 „ a„ + ^02'' "'y^i 



dt dr^ r dr 4il2 



87mV2\ip{r;t)(" 



ip{r;t)^0 



(6) 



with normalization 47r L r'^dr\(p{r;t)\'^ = 1. Here, as in 
Eq. |SJ|, d{t) is a strength parameter which is to be re- 
duced to from 1 when the radial trap is switched off. 

Equations (O and (jHJl represent nonrotating {L = 0) 
and vortex solitons {L ^ 0) rotating around z axis in the 
quasi two-dimensional and spherically-symmetric three- 
dimensional cases, respectively. To study the solitons we 
finally set d{t) = in these equations. It should be noted 
that fl is supposed to be a constant reference frequency 
and not the trap frequencies cOp or lOz. However, we took 
in Eqs. (3) — (5) the condition ujz = ^- This does not 
correspond to any specialization but only simplifies the 
equations algebraically. Nevertheless, it is of advantage 
to take 51 to have the same order of magnitude as an 
experimental trap frequency, for example 57 = 27r x 80 
Hz. With this value of 51 the dimensionless time unit 
corresponds to 57"^ = l/(27r x 80) s ~ 2 ms. In Eq. ((HJ, 
51 is a constant and fl j^ ujr- 

We solve the GP equations © and © numerically 
using the split-step time-iteration method employing 
the Crank-Nicholson discretization scheme described re- 
cently 16]. The time iteration is started with the known 
oscillator solution of these equations with zero nonlinear- 
ity n. Then in the course of time iteration the attractive 
nonlinearity is switched on very slowly and in the ini- 
tial stage the harmonic trap is also switched off slowly 
by changing d(t) from 1 to 0. If the nonlinearity is in- 
creased rapidly the system collapses. The tendency to 
collapse must be avoided to obtain a stabilized soliton. 
After switching off the harmonic trap in Eqs. ^ and © 
and after slowly introducing a final attractive nonlinear- 
ity no, if n is replaced by no[l — csin(ci;i)] a stabilization of 
the final solution could be obtained for a suitably chosen 
c and a large w. The stabilization could be obtained for 
a range of values of c and u provided that uq is negative 
corresponding to attraction. After some experimentation 
with Eqs. ^ and we opted for the choice c = 4 and 
Lu — IOtt in all our calculations except in Fig. 1 (b) in 
2D and 3D — variational and numerical. In Fig. 1 (b) 
we report some results for lo — 20n for comparison. 



III. VARIATIONAL RESULTS 
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To understand how the stabilization can take place we 
employ a variational method with the following Gaussian 



wave function for the solution of Eqs. |SJ and © [SlITOJI 



ip{r,t) =iV(i)r^exp 



2i?2(f) 



^mr' 



ia{t) 



(7) 

where N{t), R{t), /3(i), and a{t) are the normahzation, 
width, chirp, and phase of the sohton, respectively. In 
3D N{t) = [7r3/4i?3/2(^)]-i and L = 0, and in 2D 
N{t) = [7ri/2i?^+i(i)^]-i. The Lagrangian density 
for generating Eq. 101 with tj^ = is 



C{ip) = 



i (dip , dip* 



dt 



'•-P 



dt 



■ ip 



dip 



dr 



L^W? 



-^am 



(8) 

where g = 87rn\/2 in 3D and g = 4nV27r in 2D. The trial 
wave function Q is substituted in the Lagrangian density 
and the effective Lagrangian is calculated by integrating 
the Lagrangian density: Lg-g- = J C{ip)df. 

The Euler-Lagrange equations for R{t) and /3(i) are 
then obtained from the effective Lagrangian in standard 
fashion in 3D: 



dR{t) 

dt 
d(3{t) 

dt 



2R{t)P{t) 
2 



RHt) 



-2/?2(t) + 



/2^i?5(i) 



(9) 
(10) 



From Eqs. Q and (|10|l we get the following second- 
order differential equation for the evolution of the width 



d^Rjt) 
df^ 



go + gi sin(wt) 



R'(t) 



(11) 



with g — go+gi sin(ci;i), where go corresponds to the con- 
stant part of the scattering length and gi to the oscillat- 
ing part. We separate R{t) into a slowly varying part A{t) 
and a rapidly varying part Bit) by R{t) — A{t) + B{t). 
Substituting this into Eq. I|ll|l and retaining terms on 
the order of w^^ in B{t) we obtain the following equa- 
tions of motion for B{t) and A{t): 



d'^B{t) _ gism{ijjt) 
dt^ 

d'^Ajt) _ 
df^ 




2V2gi{B{t)shi{Lut)) 
7r3/M5(i) '' 



where ( ) denotes time average over rapid oscillation. 

Using the solution B{t) = -gi sm{ujt)/[V2n^uj^A^{t)], 
the equation of motion for A{t) becomes 



d^A{t) 


4 

A3 

d 

dA 


.90 gl 


(12) 
.(13) 


dt^ 


V2^A^ TT^LU^A^ 

'2 50 5? ■ 




A^ 3V2^A^ 8tt3lu^A»_ 



The quantity in the square bracket in Eq. H13|l is the 
effective potential U{A) of the equation of motion: 



U{A) 




SttS^jMS ■ 



(14) 



Small oscillations around a stable configuration are pos- 
sible when there is a minimum in this effective poten- 
tial 8]. Unfortunately, this condition does not lead to 
a simple analytical solution. However, straightforward 
numerical study reveals that this effective potential has 
a minimum for a negative go corresponding to attrac- 
tion with \go\ above a critical value. For a numerical 
calculation the quantity g is taken to be of the form 
g — go + g\ sin(ti;i) = go[l ^ 4sin(cL'i)] so that g\ = — 4(7o 
with go negative (attractive). 

In Figs. 1 (a) and (b) we plot the effective potential 
C/(A) vs. A for different go for w = IOtt and w = 207r, re- 
spectively. We find that, as the value of 50 is reduced, the 
effective potential develops a minimum which gradually 
becomes deeper and deeper. The depth of the minimum 
in the effective potential increases as lo increases. For 
w = IOtt and go ~ —100 there is no minimum in the ef- 
fective potential C/(A), whereas a minimum has appeared 
for go — —200 which becomes deeper for go = —300 and 
—500. For 10 = 207r and go = —100, a minimum has al- 
ready appeared in Fig. 1 (b). In the rest of this study we 
use the frequency to = IOtt, although its actual value has 
no consequence on the calculation so long as it is large 
corresponding to rapid oscillation. A careful examina- 
tion reveals that the threshold for the minimum in the 
present case is given by 50 ~ —116 for cj = IOtt. Hence, 
for w = IOtt stabilization is not possible for go — —100, 
and it is possible for go < —116. Also there is no upper 
limit for \go\ and stabilization is possible for an arbitrar- 
ily large \go\- As the first and the third terms on the right 
hand side (rhs) of Eq. H14(l are positive, no stabilization 
is possible for a positive go corresponding to repulsion. 
We shall verify these findings by actual numerical calcu- 
lation in the following. 

Similarly, a two-dimensional soliton of Eq. 101 leads to 
the following equation of motion for the small oscillation 
of width A{t) for a general L: 



d^Ajt) 
dt^ 



d_ 
dA 



50 



47r 



1 



gf 



7rA2 



Att^lo^A^ 



(15) 



A similar equation was obtained before for L = |3|. 
The condition for small oscillations is given by go < — 47r 

or n < — ■y/7r/2, when the first term on the rhs of Eq. 
(|15|l becomes negative allowing for the possibility of a 
minimum in the effective potential in the square bracket 
resulting in stable small oscillations. 



IV. NUMERICAL RESULTS 

With this preliminary variational study we turn to a 
full numerical investigation of Eqs. lO and ^ in 2D and 
3D. As a warm up it is worthwhile to redo the numerical 
study for a L = soliton in 2D and extend it to a L = 1 
vortex soliton before considering a L = soliton in 3D. 

There could be many ways of numerical stabilization 
of the soliton. In the course of time evolution of the GP 
equation certain initial conditions are necessary for the 
stabilization of a soliton with a specific nonlinearity |no| 
above a critical value. As one requires a large (attractive) 
nonlinearity |no| for stabilization, one needs to reduce the 
harmonic trap frequency while increasing the nonlinear- 
ity |n|. Unless the trap frequency is reduced the system 
will collapse ^3 due to attraction. In other words one 
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FIG. 1: The effective potential U{A) of Eq. ^ vs. A in 
arbitrary units for go = —100, —200, —300, and —500 (from 
upper to lower curve) for (a) lo = IOtt, and (b) u) = 207r. 



allows the system to expand and simultaneously increase 
the nonlinearity |n|. During this process the harmonic 
trap is removed, and after the final nonlinearity no is at- 
tained at time io the periodically oscillating nonlinearity 
n = no[l — 4sin{107r(i — to)}] is applied for t > tg. If the 
size of the condensate is close to the desired size, a sta- 
bilization of the condensate for a large time is obtained. 
This procedure could also be followed in an experimen- 
tal attempt to stabilize a soliton. Saito et al. used 
a qualitatively similar, but quantitatively different pro- 
cedure for stabilization. The procedure of Saito et al. 
could also be applied successfully in the present context. 
The correct implementation of the above calculational 
scheme is important for stabilization. If the (attrac- 
tive) nonlinearity after switching off the harmonic trap 
is strong for the size of the condensate the system be- 
comes highly attractive in the final stage and it eventu- 
ally collapses. If the nonlinearity after switching off the 
harmonic trap is weak for its size the system becomes 
weakly attractive in the final stage and it expands to in- 
finity. The final nonlinearity has to have an appropriate 
intermediate value, decided by trial, for eventual stabi- 
lization. In Figs. 2 (a) and (b) we provide the actual time 
variation of nonlinearity \n{t)\ as well as the strength pa- 



FIG. 2: Variation of nonlinearity parameter |n(i)| and the 
strength of radial trap d{t) in (a) Eq. l|^ and (b) Eq. @ 
in the initial stage of stabilization until the desired final non- 
linearity no is attained at time io. For t > to, the oscillating 
nonlinearity n — no[l — 4sin{107r(t — to)}] is applied. In 2D, 
no = —9.1 for L = and no = —9.7 for L = 1; in 3D, 
no = -48.9. 



rameter for radial trap d(t) employed in Eqs. |(5j and 
© for 2D and 3D, respectively. A finetuning of the fi- 
nal nonlinearity \n{t)\ was needed for stabilization over a 
large interval of time as reported in this paper. 

The results of numerical calculation based on Eq. Q 
are shown in Figs. 3 (a) and (b) for the two dimensional 
soliton (L = 0) and vortex soliton (L = 1), respectively, 
where we plot the radial part of the wave function at 
times t = 0, 50, 100, 150, 200, 250, 300 after stabilization 
is obtained. With the present value of fJ = 27r x 80 Hz, 
t = 300 corresponds to 600 ms. 

Next we turn to a numerical calculation in 3D. The 
relation between the constant g considered in the varia- 
tional calculation and nonlinearity n in the GP equation 
is g = 8nn^/2 « 35n. From the variational calculation 
presented in Fig. 1 (a) for w = IOtt, we find that the 
condition for stabilization of a soliton in 3D is go < —116 
which corresponds to no < —3.3, approximately. From a 
complete numerical solution of the GP equation © we 
also find that there is a threshold of nonlinearity for sta- 



3 0.006 

9- 



80 100 




— 0.004 



40 80 120 160 200 
P 



FIG. 3: (a) Wave function \f{p)\ of stabilized soliton in 2D 
with L = and n = — 9.1[1 — 4sin(107ri)] in Eq. ^ at times 
t = 0, 50, 100, 150, 200, 250, 300. (b) Same for vortex soliton 
with L = 1 and n = — 9.7[1 — 4sin(107ri)] in Eq. ^ at times 
t = 0, 50, 100, 150, 200, 250, 300. 



bilization consistent with the variational calculation. In 
the numerical calculation it was difficult to obtain accu- 
rately the threshold value of uq for stabilization. How- 
ever, we could not obtain stabilization of the soliton for 
no > —10. The stabilization is possible for stronger no 
(larger |no|). In Fig. 4 we plot the radial wave function of 
Eq. ® in 3D at times t = 0, 50, 100, 150, 200, 250, 300 for 
no = —48.9 after obtaining the stabilization. The nar- 
row spread of the wave function over the large interval of 
time shows the quality of stabilization. In both Figs. 3 
and 4 the results at intermediate times lie in the region 
covered by the plots. The plot of the full wave function, 
rather than that of mean radii or the wave function at 
a particular point vs time, clearly shows the degree of 
stabilization achieved. 

Of the three cases presented in Figs. 3 and 4 the vortex 
soliton of Fig. 3 (b) is the most stable with minimum 
oscillation. For the rotating (vortex) soliton, the outward 
centrifugal force approximately balances the attractive 
inward force (as seen in a rotating frame). The final 
exact balance is provided by the oscillating nonlinearity. 
However, for L = there is no centrifugal force and the 
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FIG. 4: Wave function [(/^(r)! of stabilized soliton in 3D 
with n = -48.9[1 - 4sin(107ri)] in Eq. © at times t = 
0,50,100,150,200,250,300. 



stabilization is more difhcult. 

The stabilization in both 2D and 3D can be obtained 
for (attractive) BEG solitons with nonlinearity |no| larger 
than a critical value. In 2D the variational critical non- 
linearity is n^j.^^ = — ■y/7r/2, whereas the final average 
nonlinearities in Figs. 3 (a) and (b) are no = —9.1 and 
—9.7, respectively. In actual numerical calculation we 
found that stronger the nonlinearity |no|, more sustained 
was the stabilization of the soliton. The effective poten- 
tial develops a deeper minimum for a larger nonlinearity 
no|. The variational threshold for stabilization in 3D is 
'^critl ~ 3.3. In Fig. 4 the actual final average non- 
linearity I no I = 48.9 is much larger than the variational 



threshold 



'crit I 



3.3. 



Using a variational procedure alone, not quite iden- 
tical with the present approach, AbduUaev et al. [l(l | 
also had found that a stabilization of a soliton could be 
possible in 3D via a temporal modulation of the scat- 
tering length. However, they confirmed after further an- 
alytical and numerical studies that such a stabilization 
does not take place in 3D. Saito et al. Q, on the other 
hand, are silent about the possibility of the stabiliza- 
tion of a soliton in 3D. We point out one possible reason 
for the negative result obtained in 3D [l3|. The non- 
linearity parameter A = — 1 used in Ref. |To| for sta- 
bilizing a soliton in 3D corresponds in our notation to 
no = -V2^/[87rV2] = -\/^8 « -0.22. [The relation 
between no and A of Ref. [lO^ follows from present Eq. 
(HU and their Eq. (37)]. The nonhnearity A = -1 is 
much too weak for obtaining a stabilized soliton in 3D. 
It should be noted that the value of no used for the sta- 
bilization of a soliton in 3D in the present calculation is 
—48.9, whereas the variational threshold for stabilization 
is nj,pj^ w —3.3. These values of nonlinearities are much 
stronger than the value no = —0.22 used in Ref. |10| . 



V. DISCUSSION AND CONCLUSION 

In this paper we have discussed the stabihzation of 
a bright vortex sohton in 2D and a bright soUton in 
3D by a periodic temporal modulation of the scattering 
length. Now we compare the dimensionless parameters 
used in the simulations to typical numbers for an exper- 
imental system ^^Rb. This system has a Feshbach reso- 
nance which can be used to vary the effective interaction 
betw een atoms by varying the atomic scattering length 
\n\. With the reference trap f] = 27r x 80 Hz, for ^^^Rb 

the dimensionless length parameter I = y/K/JmTf) w 1.2 
/xm. In a 3D condensate of 10000 atoms, the variational 
critical nonlinearity l^icritl ^ "^"^ ^^^ ^^^^ system corre- 
sponds to a scattering length a = ri^^^-^l/N « —0.4 nm. 
The applied nonlinearity ng = —48.9 of the present 3D 
simulation corresponds to a scattering length a « —5.9 
nm. The oscillating nonlinearity corresponds in this case 
(a w —5.9 nm) to a variation of scattering length be- 
tween —30 nm and 18 nm. For 100000 atoms the above 
values for scattering length will be reduced by a factor of 
10. Similar variations of scattering length of ^^Rb have 
already been realized in the laboratory via a Feshbach 
resonance in actual experiments 17]. Hence it might be 
possible to stabilize a ^^Rb condensate using a Feshbach 
resonance. With fi = 27r x 80 Hz, the interval of stabi- 
lization of 300 units of time in Figs. 3 and 4 corresponds 
to 600 ms, which is a reasonably large interval of time. 

The present study has important consequence in the 
generation of a stable spatiotemporal soliton of nonlin- 
ear optics 0, 18] which is an optical wave packet con- 
fined in all three directions and often referred to as a 
light bullet. The existence and stability of self trapped 
beams in a nonlinear medium is a subject of active re- 



search since its suggestion ^3- Such a spatiotemporal 
optical soliton satisfies an equation in an anomalously 
dispersive medium quite similar to Eq. ^ with ujr — 
^3]. Hence a stable solution for a light bullet in actual 
3D can be obtained through a modulation of the cubic 
Kerr nonlinearity. This modulation can be achieved in a 
layered nonlinear medium with sign-altering Kerr nonlin- 
earity 10, 13]. The possibility of such a stabilization in 
two space dimensions was demonstrated in 13], whereas 
its impossibility in 3D has been emphasized in jTj ] . 

In conclusion, from a numerical solution of the GP 
equation we find that it is possible to stabilize a matter- 
wave bright soliton in 3D and a vortex soliton in 2D by 
employing a rapid periodic modulation of the scatter- 
ing length a via a Feshbach resonance with an attractive 
(negative) mean value ao via: a —^ ao[l — cshi(ujt)] with 
a large c and uj. From a variational calculation we show 
that this oscillation produces a minimum in the effec- 
tive potential, thus producing a potential well in which 
the soliton can be trapped and execute small oscillations. 
The sinusoidal variation of a is actually not needed for 
stabilization, any periodic fluctuation between positive 
and negative values stabilizes the soliton. This is of in- 
terest to investigate if such BEG "bullets" could be cre- 
ated experimentally in 3D. As the mathematical equa- 
tion satisfied by a light bullet 2] in 3D is quite similar to 
the nonlinear three-dimensional equation © we suggest 
the possibility of creating light bullets in a layered Kerr 
medium with sign-altering nonlinearity. 
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